// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtUintahFileFormat.h                          //
// ************************************************************************* //

#include <avtUintahFileFormat.h>
#include <avtUintahOptions.h>

#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkCellData.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkLongArray.h>
#include <vtkIntArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtParallel.h>
#include <avtDatabase.h>
#include <avtDatabaseMetaData.h>
#include <avtIntervalTree.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtStructuredDomainNesting.h>
#include <avtVariableCache.h>
#include <avtMaterial.h>
#include <avtCallback.h>

#include <TimingsManager.h>
#include <DebugStream.h>
#include <DBOptionsAttributes.h>
#include <InvalidDBTypeException.h>
#include <InvalidFilesException.h>
#include <InvalidVariableException.h>
#include <InstallationFunctions.h>

#include <dlfcn.h>
#include <string>
#include <vector>

#ifdef PARALLEL
#include <mpi.h>
// Only get the option for serialized reads if compiling the parallel version
//#define SERIALIZED_READS
#endif


using namespace UintahDBOptions;

const double NAN_REPLACE_VAL = 1.0E9;

// ****************************************************************************
//  Method: avtUintahFileFormat Constructor
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
avtUintahFileFormat::avtUintahFileFormat(const char *filename,
                                         const DBOptionsAttributes* attrs) :
  avtMTMDFileFormat(filename)
{
  // int t1 = visitTimer->StartTimer();

  for (int i=0; attrs!=0 && i<attrs->GetNumberOfOptions(); ++i)
  {
    if (attrs->GetName(i) == UINTAH_LOAD_EXTRA)
    {
      loadExtraGeometry = attrs->GetEnum(UINTAH_LOAD_EXTRA);
    }
    else if (attrs->GetName(i) == UINTAH_DATA_VARIES_OVER_TIME)
    {
      dataVariesOverTime = attrs->GetBool(UINTAH_DATA_VARIES_OVER_TIME);
    }
  }

  // Verify that it is a UDA index.xml file:
  // The 2nd line should look like this <Uintah_DataArchive>.
  FILE * fp = fopen( filename, "r" );
  if( fp == NULL )
  {
    std::string error = std::string( "Failed to open file: " ) + filename;
    EXCEPTION1( InvalidDBTypeException, error.c_str() );
  }

  char line[1024];
  char * result = fgets( line, 1024, fp );
  if( result )
  {
    result = fgets( line, 1024, fp );
  }

  std::string lineStr = line;
  if( !result || lineStr.find( "<Uintah_DataArchive>" ) == std::string::npos )
  {
    std::string error = std::string( filename ) +
      " does not appear to be a <Uintah_DataArchive>.";
//    printf("here: %s\n", error.c_str());
    EXCEPTION1( InvalidDBTypeException, error.c_str() );
  }

  fclose( fp );

  // This environment variable prevents Core/Thread/Thread_pthreads.cc
  // from adding it's atexit handler, which will kill visit's process
  // when dlclose is called.
  char *tmp = new char[50];
  strcpy(tmp, "THREAD_NO_ATEXIT=1");
  putenv(tmp);

  int dlopen_mode = RTLD_NOW;

#ifdef UINTAH_INTERFACES_LIB
  const char * lib_name = UINTAH_INTERFACES_LIB;
#else
#error  "UINTAH_INTERFACES_LIB has not been defined"
  const char * lib_name = "UINTAH_INTERFACES_LIB";
#endif

  setenv("THREAD_NO_CATCH_SIGNALS", "True", 1);

  // First try to open the library without any paths - assumes the
  // user has the Uintah library parth in their *_LIBRARY_PATH.
  libHandle = dlopen(lib_name, dlopen_mode);

  // if( libHandle )
  //   std::cerr << __LINE__ << " Uintah lib " << lib_name << std::endl;

  // Try the visit installation library directory.
  if (!libHandle)
  {
    const char *vlibdir = GetVisItLibraryDirectory().c_str();

    if( vlibdir )
    {
      char *lib = (char *) malloc( strlen(vlibdir) + strlen(lib_name) + 8 );

      sprintf( lib, "%s/%s", vlibdir, lib_name );

      libHandle = dlopen(lib, dlopen_mode);

      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
    }

    // Try a relative installed path
    if (!libHandle)
    {
      char *lib = (char *) malloc( strlen(vlibdir) + strlen(lib_name) + 12 );

      sprintf( lib, "%s/uintah/%s", vlibdir, lib_name );

      libHandle = dlopen(lib, dlopen_mode);

      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
    }
  }

  // Did not open so try to open the library using the path from the
  // calling library which would be in plugins/databases.
#ifdef HAVE_WINDOWS_H
  // Not sure to do here???
#else
  if (!libHandle)
  {
    char *pathname = 0;

    Dl_info info;
    if (dladdr(__builtin_return_address(0), &info))
    {
      const char *lastslash = strrchr(info.dli_fname,'/');

      if( lastslash )
      {
        int pathLen = strlen(info.dli_fname) - strlen(lastslash);

        pathname = (char *) malloc(pathLen+2);
        strncpy( pathname, info.dli_fname, pathLen );
        pathname[pathLen] = '\0';
      }
    }

    if( pathname )
    {
      char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 8 );

      sprintf( lib, "%s/%s", pathname, lib_name );

      libHandle = dlopen(lib, dlopen_mode);

      // if( libHandle )
      //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;

      // Try a relative installed path
      if (!libHandle)
      {
        char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 24 );

        sprintf( lib, "%s/../../lib/uintah/%s", pathname, lib_name );

        libHandle = dlopen(lib, dlopen_mode);

        // if( libHandle )
        //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
      }

      // Try a relative installed path
      if (!libHandle)
      {
        char *lib = (char *) malloc( strlen(pathname) + strlen(lib_name) + 16 );

        sprintf( lib, "%s/../../lib/%s", pathname, lib_name );

        libHandle = dlopen(lib, dlopen_mode);

        // if( libHandle )
        //        std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
      }

      free( pathname );
    }
  }

#endif

  // Did not open so try to open the library using the library path
  // found where the user installed Uintah.
#if __APPLE__
  if (!libHandle)
  {
    char *lib = (char *) malloc( strlen(lib_name) + 32 );

    sprintf( lib, "@executable_path/../lib/uintah/%s", lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

  if (!libHandle)
  {
    char *lib = (char *) malloc( strlen(lib_name) + 32 );

    sprintf( lib, "@executable_path/../lib/%s", lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

#endif

  // Did not open so try to open the library using the Uintah library path
  // found when VisIt was configured.
  if (!libHandle)
  {
#ifdef UINTAH_LIBRARY_DIR
    const char * lib_path = UINTAH_LIBRARY_DIR;
#else
#error  "UINTAH_LIBRARY_DIR has not been defined"
#endif

    char *lib = (char *) malloc( strlen(lib_path) + strlen(lib_name) + 8 );

    sprintf( lib, "%s/%s", lib_path, lib_name );

    libHandle = dlopen(lib, dlopen_mode);

    // if( libHandle )
    //   std::cerr << __LINE__ << " Uintah lib " << lib << std::endl;
  }

  // Library was never opened so report back an error.
  if (!libHandle)
  {
    char* errString = dlerror();
    std::string str = "The library " + std::string(lib_name) +
      " could not be opened. The following error was reported: " +
      std::string( errString );
    EXCEPTION1(InvalidDBTypeException, str.c_str());
  }

  // All possible function calls - check here
  char *error;

  openDataArchive = (DataArchive* (*)(const std::string&)) dlsym(libHandle, "openDataArchive");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function openDataArchive could not be located in the library!!!");
  }

  closeDataArchive = (void (*)(DataArchive*)) dlsym(libHandle, "closeDataArchive");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function closeDataArchive could not be located in the library!!!");
  }

  getGrid = (GridP* (*)(DataArchive*, int)) dlsym(libHandle, "getGrid");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getGrid could not be located in the library!!!");
  }

  releaseGrid = (void (*)(GridP*)) dlsym(libHandle, "releaseGrid");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function releaseGrid could not be located in the library!!!");
  }

  queryProcessors = (unsigned int (*)(DataArchive*)) dlsym(libHandle, "queryProcessors");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function queryProcessors could not be located in the library!!!");
  }

  getCycleTimes = (std::vector<double> (*)(DataArchive*)) dlsym(libHandle, "getCycleTimes");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getCycleTimes could not be located in the library!!!");
  }

  getTimeStepInfo = (TimeStepInfo* (*)(DataArchive*, GridP*, int, int)) dlsym(libHandle, "getTimeStepInfo");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getTimeStepInfo could not be located in the library!!!");
  }

  getGridData = (GridDataRaw* (*)(DataArchive*, GridP*, int, int, std::string, int, int, int[3], int[3], int)) dlsym(libHandle, "getGridData");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getGridData could not be located in the library!!!");
  }

#if LIB_VERSION_LE(UINTAH, 2,0,0)
  variableExists = (bool (*)(DataArchive*, std::string)) dlsym(libHandle, "variableExists");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function variableExists could not be located in the library!!!");
  }
#endif

#if LIB_VERSION_LE(UINTAH, 2,5,1)
  getNumberParticles = (unsigned int (*)(DataArchive*, GridP*, int, int, int, int)) dlsym(libHandle, "getNumberParticles");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getNumberParticles could not be located in the library!!!");
  }
#endif

  getParticleData = (ParticleDataRaw* (*)(DataArchive*, GridP*, int, int, std::string, int, int)) dlsym(libHandle, "getParticleData");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getParticleData could not be located in the library!!!");
  }

  getParticlePositionName = (std::string (*)(DataArchive*)) dlsym(libHandle, "getParticlePositionName");
  if((error = dlerror()) != NULL)
  {
    EXCEPTION1(InvalidDBTypeException, "The function getParticlePositionName could not be located in the library!!!");
  }

  // Use the folder name, not the index.xml file name to open the archive
  std::string folder(filename);
  size_t found = folder.find_last_of("/");
  folder = folder.substr(0, found);
  archive = (*openDataArchive)(folder);

  // Timestep times
  // int t2 = visitTimer->StartTimer();

  nProcs = (*queryProcessors)(archive);
  cycleTimes = (*getCycleTimes)(archive);

  // visitTimer->StopTimer(t2, "avtUintahFileFormat::avtUintahFileFormat getCycleTimes");

  // Haven't loaded any timestep data yet
  stepInfo = NULL;
  currTimeStep = -1;

  // visitTimer->StopTimer(t1, "avtUintahFileFormat::avtUintahFileFormat");
}


// ****************************************************************************
//  Method: avtUintahFileFormat Destructor
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************

avtUintahFileFormat::~avtUintahFileFormat()
{
  if (grid)
    (*releaseGrid)(grid);

  if (archive)
    (*closeDataArchive)(archive);

  if (stepInfo)
    delete stepInfo;

  dlclose(libHandle);
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
int
avtUintahFileFormat::GetNTimesteps(void)
{
  return cycleTimes.size();
}


// ****************************************************************************
// Method: avtUintahFileForma::GetTime
//
// Purpose:
//   Get the time.
//
// Programmer: sshankar
// Creation:   Fri Feb 6 15:31 MST 2009
//
// ****************************************************************************
double
avtUintahFileFormat::GetTime(int ts)
{
  return cycleTimes[ts];
}


// ****************************************************************************
// Method: avtUintahFileForma::ActivateTimestep
//
// Purpose:
//   Get ready to read data for the given timestep
//
// ****************************************************************************
void
avtUintahFileFormat::ActivateTimestep(int ts)
{
  if (currTimeStep == ts)
    return;

  // int t1 = visitTimer->StartTimer();

  // Get the uda grid for the new timestep
  // int t2 = visitTimer->StartTimer();

  if (grid)
    (*releaseGrid)(grid);

  grid = (*getGrid)(archive, ts);

  // visitTimer->StopTimer(t2, "avtUintahFileFormat::ActivateTimestep getGrid");

  // Get the time step info for the new timestep
  // int t3 = visitTimer->StartTimer();

  if (stepInfo)
    delete stepInfo;

  stepInfo = (*getTimeStepInfo)(archive, grid, ts, loadExtraGeometry);

  currTimeStep = ts;
  forceMeshReload = true;

  // visitTimer->StopTimer(t3, "avtUintahFileFormat::ActivateTimestep getTimeStepInfo");

  // visitTimer->StopTimer(t1, "avtUintahFileFormat::ActivateTimestep");
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::AddExpressionsToMetadata
//
//  Purpose:
//      Adds the Uintah specific expressions to the meta data.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
void
avtUintahFileFormat::AddExpressionsToMetadata(avtDatabaseMetaData *md)
{
  std::string home(getenv("HOME"));
  std::string fname = home +"/.visit/udaExpressions.txt";

  FILE *f = fopen(fname.c_str(), "r");
  if (!f)
  {
    debug5<<"couldn't open uda expressions file"<<endl;
    return;
  }

  int line = 0;
  char cline[2048];
  while (fgets(cline, 2048, f))
  {
    line++;

    std::string line(cline);
    if (line[0] == '#' ||
        line.length() == 0 ||
        line == "\n" || line=="\r\n")
      continue;

    // first token is the expression name
    int space = line.find(" ");
    if (space < 0)
    {
      debug5<<"uda expression syntax error on line "<<line<<endl;
      continue;
    }

    std::string name = line.substr(0, space);
    line = line.substr(space+1);

    // second token is the expression type
    space = line.find(" ");
    if (space < 0)
    {
      debug5<<"uda expression syntax error on line "<<line<<endl;
      continue;
    }

    std::string type = line.substr(0, space);
    line = line.substr(space+1);

    Expression::ExprType etype = Expression::Unknown;
    if (type == "scalar")
      etype = Expression::ScalarMeshVar;
    if (type == "vector")
      etype = Expression::VectorMeshVar;
    if (type == "tensor")
      etype = Expression::TensorMeshVar;
    if (type == "symtensor")
      etype = Expression::SymmetricTensorMeshVar;
    if (type == "array")
      etype = Expression::ArrayMeshVar;
    if (type == "curve")
      etype = Expression::CurveMeshVar;

    if (etype == Expression::Unknown)
    {
      debug5<<"uda expression - unknown type on line "<<line<<endl;
      continue;
    }


    Expression e;
    e.SetName(name);
    e.SetType(etype);
    e.SetDefinition(line);
    md->AddExpression(&e);
  }

  fclose(f);
}


// ****************************************************************************
//  Method: avtUintahFileFormat::ReadMetaData
//
//  Purpose:
//      Does the actual work for PopulateMetaData().
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
void
avtUintahFileFormat::ReadMetaData(avtDatabaseMetaData *md, int timeState)
{
  // int t1 = visitTimer->StartTimer();

  ActivateTimestep(timeState);

  // Don't add node data unless NC_Mesh exists (some only have CC_MESH
  // or SFC*_MESH).
  bool addNodeData = false;

  // Don't add patch data unless CC_Mesh or NC_Mesh exists (some only
  // have SFC*_MESH).
  bool addPatchData = false;

  // Grid meshes are shared between materials, and particle meshes are
  // shared between variables - keep track of what has been added so
  // they're only added once.
  std::set<std::string> meshes_added;

  // If a variable exists with multiple materials, don't add it more
  // than once to the meta data, as it can mess up VisIt's expressions
  // variable lists.
  std::set<std::string> mesh_vars_added;

  // Loop through all vars and add them to the meta data.
  for (int i=0; i<(int)stepInfo->varInfo.size(); ++i)
  {
    // Nothing needs to be modifed for particle data, as they exist
    // only on a single level
    if (stepInfo->varInfo[i].type.find("ParticleVariable") != std::string::npos)
    {
      std::string varname = stepInfo->varInfo[i].name;
      std::string vartype = stepInfo->varInfo[i].type;

      if (vartype.find("filePointer") != std::string::npos)
      {
        continue;
      }

      // j=-1 -> all materials (*)
      for (int j=-1; j<(int)stepInfo->varInfo[i].materials.size(); ++j)
      {
        std::string mesh_for_this_var = std::string("Particle_Mesh/");
        std::string newVarname = varname+"/";

        if (j >= 0)
        {
          char buffer[128];
          sprintf(buffer, "%d", stepInfo->varInfo[i].materials[j]);
          mesh_for_this_var.append(buffer);
          newVarname.append(buffer);
        }
        else
        {
          mesh_for_this_var.append("*");
          newVarname.append("*");
        }

        addParticleMesh( md, meshes_added, mesh_for_this_var, stepInfo );

        addMeshVariable( md, mesh_vars_added,
                         newVarname, vartype, mesh_for_this_var, AVT_NODECENT );
      }
    }
    // Volume Data.
    else
    {
      bool isPerPatchVar = false;

      std::string varname = stepInfo->varInfo[i].name;
      std::string vartype = stepInfo->varInfo[i].type;

      avtCentering cent;
      std::string mesh_for_this_var;

      if (vartype.find("NC") != std::string::npos)
      {
        mesh_for_this_var.assign("NC_Mesh");
        cent = AVT_NODECENT;

        addNodeData = true;
        addPatchData = true;
      }
      else if (vartype.find("CC") != std::string::npos)
      {
        mesh_for_this_var.assign("CC_Mesh");
        cent = AVT_ZONECENT;

        addPatchData = true;
      }
      else if (vartype.find("SFC") != std::string::npos)
      {
        if (vartype.find("SFCX") != std::string::npos)
          mesh_for_this_var.assign("SFCX_Mesh");
        else if (vartype.find("SFCY") != std::string::npos)
          mesh_for_this_var.assign("SFCY_Mesh");
        else if (vartype.find("SFCZ") != std::string::npos)
          mesh_for_this_var.assign("SFCZ_Mesh");

        cent = AVT_ZONECENT;
      }
      else if (vartype.find("PerPatch") != std::string::npos)
      {
        if (varname.find("FileInfo") == 0 ||
            varname.find("CellInformation") == 0 ||
            varname.find("CutCellInfo") == 0)
          continue;

        mesh_for_this_var = "Patch_Mesh";
        cent = AVT_ZONECENT;

        isPerPatchVar = true;
      }
      else if (vartype.find("ReductionVariable") != std::string::npos ||
               vartype.find("SoleVariable")      != std::string::npos)
      {
        continue;
      }
      else
      {
        std::stringstream msg;
        msg << "Visit - "
            << "Uintah variable \"" << varname << "\"  "
            << "has an unknown variable type \""
            << vartype << "\"";

        avtCallback::IssueWarning(msg.str().c_str());
        debug5 << msg.str() << std::endl;

        continue;
      }

      addRectilinearMesh( md, meshes_added, mesh_for_this_var, stepInfo );

      // Add mesh vars for each material. If no materials loop just once.
      int numMaterials = stepInfo->varInfo[i].materials.size();

      for (int j=0; j<std::max(1,numMaterials); ++j)
      {
        std::string newVarname = varname;

        // PerPatch vars do not have a material index.
        if( isPerPatchVar )
        {
          newVarname = "Patch/" + newVarname;
        }
        // Add the material index.
        else if( numMaterials )
          {
            char buffer[128];
            sprintf(buffer, "%d", stepInfo->varInfo[i].materials[j]);
            newVarname.append("/");
            newVarname.append(buffer);
          }
        // For variables with no materials add a material index of 0.
        else
          newVarname.append("/0");

        addMeshVariable( md, mesh_vars_added,
                         newVarname, vartype, mesh_for_this_var, cent );
      }
    }
  }

  // Add the node data
  if (addNodeData)
  {
    avtVectorMetaData *vector = new avtVectorMetaData();
    vector->name = "Patch/Nodes";
    vector->meshName = "NC_Mesh";
    vector->centering = AVT_NODECENT;
    vector->hasDataExtents = false;
    vector->varDim = 3;
    md->Add(vector);
  }

  // Add the patch data
  if (addPatchData)
  {
    std::string mesh_for_this_var = "Patch_Mesh";

    addRectilinearMesh( md, meshes_added, mesh_for_this_var, stepInfo );

    avtCentering cent = AVT_ZONECENT;

    avtScalarMetaData *scalar;

#if LIB_VERSION_LE(UINTAH, 2,1,0)
    scalar = new avtScalarMetaData();
    scalar->name = "Patch/Id";
    scalar->meshName = mesh_for_this_var;
    scalar->centering = cent;
    scalar->hasDataExtents = false;
    scalar->treatAsASCII = false;
    md->Add(scalar);
#endif

    scalar = new avtScalarMetaData();
    scalar->name = "Patch/ProcId";
    scalar->meshName = mesh_for_this_var;
    scalar->centering = cent;
    scalar->hasDataExtents = false;
    scalar->treatAsASCII = false;
    md->Add(scalar);

    for (std::set<std::string>::iterator it=meshes_added.begin();
         it!=meshes_added.end(); ++it)
    {
      if ( (*it).find("NC") != std::string::npos ||
           (*it).find("CC") != std::string::npos )
      {
        avtVectorMetaData *vector = new avtVectorMetaData();
        vector->name = "Patch/Bounds/Low/" + *it;
        vector->meshName = mesh_for_this_var;
        vector->centering = cent;
        vector->hasDataExtents = false;
        vector->varDim = 3;
        md->Add(vector);

        vector = new avtVectorMetaData();
        vector->name = "Patch/Bounds/High/" + *it;
        vector->meshName = mesh_for_this_var;
        vector->centering = cent;
        vector->hasDataExtents = false;
        vector->varDim = 3;
        md->Add(vector);
      }
    }
  }

  int numLevels = stepInfo->levelInfo.size();

  int totalPatches = 0;
  for (int i = 0; i < numLevels; i++)
    totalPatches += stepInfo->levelInfo[i].patchInfo.size();

  // Used for the domain partitioning
  std::vector<int> groupIds(totalPatches);

  for (int i = 0; i < totalPatches; i++)
  {
    int level, local_patch;

    GetLevelAndLocalPatchNumber(i, level, local_patch);
    groupIds[i] = level;
  }

  md->AddGroupInformation(numLevels, totalPatches, groupIds);
  md->AddDefaultSILRestrictionDescription(std::string("!TurnOnAll"));

  // Set the cycles and times.
  md->SetCyclesAreAccurate(true);

  std::vector<int> cycles;

  cycles.resize( cycleTimes.size() );

  for(int i=0; i<(int)cycleTimes.size(); ++i )
    cycles[i] = i;

  md->SetCycles( cycles );

  md->SetTimesAreAccurate(true);
  md->SetTimes( cycleTimes );

  AddExpressionsToMetadata(md);

  // visitTimer->StopTimer(t1, "avtUintahFileFormat::ReadMetaData");
}


// ****************************************************************************
//  Method: avtUintahFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
void
avtUintahFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                              int timeState)
{
#ifdef SERIALIZED_READS
  int numProcs, rank;
  int msg = 128, tag = 256;
  MPI_Status status;

  MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
  MPI_Comm_rank(VISIT_MPI_COMM, &rank);
  //debug5 << "Proc: " << rank << " sent to mdserver" << endl;

  if (rank == 0)
  {
    ReadMetaData(md, timeState);
    MPI_Send(&msg, 1, MPI_INT, 1, tag, VISIT_MPI_COMM);
  }
  else
  {
    MPI_Recv(&msg, 1, MPI_INT, rank - 1, tag, VISIT_MPI_COMM, &status);
    if (msg == 128 && tag == 256)
    {
      ReadMetaData(md, timeState);
      if (rank < (numProcs - 1))
        MPI_Send(&msg, 1, MPI_INT, rank + 1, tag, VISIT_MPI_COMM);
    }
  }
#else
  ReadMetaData(md, timeState);
#endif
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetLevelAndLocalPatchNumber
//
//  Purpose:
//      Translates the global patch identifier to a refinement level and patch
//      number local to that refinement level.
//
//  Programmer: sshankar, taken from implementation of the plugin, CHOMBO
//  Creation:   May 20, 2008
//
// ****************************************************************************
void
avtUintahFileFormat::GetLevelAndLocalPatchNumber(int global_patch,
                                                 int &level, int &local_patch)
{
//  int t1 = visitTimer->StartTimer();

  int num_levels = stepInfo->levelInfo.size();
  int num_patches = 0;
  local_patch = global_patch;
  level = 0;

  while (level < num_levels)
  {
    num_patches = stepInfo->levelInfo[level].patchInfo.size();

    if (local_patch < num_patches)
      break;

    local_patch -= num_patches;
    level++;
  }

  //  visitTimer->StopTimer(t1, "avtUintahFileFormat::GetLevelAndLocalPatchNumber");
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetGlobalDomainNumber
//
//  Purpose:
//      Translates the level and local patch number into a global patch id.
//
// ****************************************************************************
int
avtUintahFileFormat::GetGlobalDomainNumber(int level, int local_patch)
{
  int global_patch = 0;

  for (int l=0; l<level; ++l)
    global_patch += stepInfo->levelInfo[l].patchInfo.size();

  global_patch += local_patch;

  return global_patch;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetDomainBoundariesAndNesting
//
//  Purpose:
//      Calculates two important data structures.  One is the structure domain
//      nesting, which tells VisIt how the AMR patches are nested, which allows
//      VisIt to ghost out coarse zones that are refined by smaller zones.
//      The other structure is the rectilinear domain boundaries, which tells
//      VisIt which patches are next to each other, allowing VisIt to create
//      a layer of ghost zones around each patch.  Note that this only works
//      within a refinement level, not across refinement levels.
//
//
// NOTE: The cache variable for the mesh MUST be called "any_mesh",
// which is a problem when there are multiple meshes or one of them is
// actually named "any_mesh" (see https://github.com/visit-dav/visit/issues/138
// Thus, for each mesh we keep around our own cache variable and if this
// function finds it then it just uses it again instead of recomputing it.
//
// ****************************************************************************
void
avtUintahFileFormat::GetDomainBoundariesAndNesting(int timestate,
                                                   const std::string &meshname)
{
#ifdef MDSERVER
    return;

#else
  // int t1 = visitTimer->StartTimer();

  // Lookup the mesh in the cache and if it's not there, compute it
  if (*this->mesh_domains[meshname] == NULL || forceMeshReload == true)
  {
    // Calculate some info needed in the rest of the routine.
    int num_levels = stepInfo->levelInfo.size();
    int totalPatches = 0;
    for (int level = 0 ; level < num_levels ; level++)
      totalPatches += stepInfo->levelInfo[level].patchInfo.size();

    // Set up the data structure for the patch boundaries. The data
    // does all the work. It just needs to know the extents of each
    // patch.
    avtRectilinearDomainBoundaries *rdb =
      new avtRectilinearDomainBoundaries(true);

    rdb->SetNumDomains(totalPatches);

    //debug5<<"Calculating avtRectilinearDomainBoundaries for "<<meshname<<" mesh ("<<rdb<<").\n";

    for (int patch=0; patch<totalPatches; patch++)
    {
      int my_level, local_patch;
      GetLevelAndLocalPatchNumber(patch, my_level, local_patch);

      PatchInfo &patchInfo =
        stepInfo->levelInfo[my_level].patchInfo[local_patch];

      int plow[3], phigh[3];
      patchInfo.getBounds(plow, phigh, meshname);

      // For node based meshes add one if there is a neighbor patch.
      if( meshname.find("NC_") == 0 )
      {
        int nlow[3], nhigh[3];
        patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

        for (int i=0; i<3; i++)
          phigh[i] += nhigh[i];
      }

      // These are indices, the high values are exclusive.
      int extents[6] = { plow[0], phigh[0],
                         plow[1], phigh[1],
                         plow[2], phigh[2] };

      // These extents are exclusive.
      rdb->SetIndicesForAMRPatch(patch, my_level, extents);

      // debug5 << "\trdb->SetIndicesForAMRPatch(" << patch << "," << my_level
      //        << ", <" << extents[0] << "," << extents[2] << "," << extents[4]
      //        << "> to <" << extents[1] << "," << extents[3] << "," << extents[5]
      //        << ">)\n";
    }

    rdb->CalculateBoundaries();

    this->mesh_boundaries[meshname] =
      void_ref_ptr(rdb, avtStructuredDomainBoundaries::Destruct);


    // Domain Nesting
    avtStructuredDomainNesting *dn =
      new avtStructuredDomainNesting(totalPatches, num_levels);
    dn->SetNumDimensions(3);

    // debug5 << "Calculating avtStructuredDomainNesting for "
    //     << meshname << " mesh (" << dn << ")." << std::endl;

    // Calculate the refinement ratio from one level to the next.
    for (int level = 0; level<num_levels; ++level)
    {
      // SetLevelRefinementRatios requires data as a vector<int>
      std::vector<int> rr(3);
      for (int i=0; i<3; i++)
        rr[i] = stepInfo->levelInfo[level].refinementRatio[i];

      dn->SetLevelRefinementRatios(level, rr);

      // debug5 << "\tdn->SetLevelRefinementRatios(" << level << ", <"
      //             << rr[0] << "," << rr[1] << "," << rr[2] << ">)\n";
    }

    // Calculating the child patches really needs some better sorting
    // than what is crrently being done.  This is likely to become a
    // bottleneck in extreme cases.  Although this routine has
    // performed well for a previous 55K patch run.
    std::vector< std::vector<int> > childPatches(totalPatches);

    for (int level=num_levels-1; level>0; --level)
    {
      int prev_level = level-1;
      LevelInfo &levelInfoParent = stepInfo->levelInfo[prev_level];
      LevelInfo &levelInfoChild = stepInfo->levelInfo[level];

      for (int child=0; child<(int)levelInfoChild.patchInfo.size(); ++child)
      {
        PatchInfo &childPatchInfo = levelInfoChild.patchInfo[child];

        int child_low[3], child_high[3];
        childPatchInfo.getBounds(child_low, child_high, meshname);

        // For node based meshes add one if there is a neighbor patch.
        if( meshname.find("NC_") == 0 )
        {
          int nlow[3], nhigh[3];
          childPatchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

          for (int i=0; i<3; i++)
            child_high[i] += nhigh[i];
        }

        int nParents = levelInfoParent.patchInfo.size();
        for (int parent = 0; parent<nParents; ++parent)
        {
          PatchInfo &parentPatchInfo = levelInfoParent.patchInfo[parent];

          int parent_low[3], parent_high[3];
          parentPatchInfo.getBounds(parent_low, parent_high, meshname);

          // For node based meshes add one if there is a neighbor patch.
          if( meshname.find("NC_") == 0 )
          {
            int nlow[3], nhigh[3];
            parentPatchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

            for (int i=0; i<3; i++)
              parent_high[i] += nhigh[i];
          }

          int mins[3], maxs[3];
          for (int i=0; i<3; ++i)
          {
            mins[i] = std::max( child_low[i],
                                parent_low[i] *levelInfoChild.refinementRatio[i]);
            maxs[i] = std::min( child_high[i],
                                parent_high[i]*levelInfoChild.refinementRatio[i]);
          }

          bool overlap = (mins[0] < maxs[0] &&
                          mins[1] < maxs[1] &&
                          mins[2] < maxs[2]);

          if (overlap)
          {
            int child_gpatch = GetGlobalDomainNumber(level, child);
            int parent_gpatch = GetGlobalDomainNumber(prev_level, parent);
            childPatches[parent_gpatch].push_back(child_gpatch);
          }
        }
      }
    }

    // Now that the extents for each patch is known and what its
    // children are, pass the structured domain boundary that
    // information.
    for (int p=0; p<totalPatches; ++p)
    {
      int my_level, local_patch;
      GetLevelAndLocalPatchNumber(p, my_level, local_patch);

      PatchInfo &patchInfo =
        stepInfo->levelInfo[my_level].patchInfo[local_patch];

      int plow[3], phigh[3];
      patchInfo.getBounds(plow, phigh, meshname);

      std::vector<int> extents(6);

      // For node based meshes add one if there is a neighbor patch.
      if( meshname.find("NC_") == 0 )
      {
        int nlow[3], nhigh[3];
        patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

        // For node meshes always subtract two because the domain is
        // inclusive.
        for (int i=0; i<3; i++)
        {
          phigh[i] += nhigh[i];

          extents[i+0] = plow[i];
          extents[i+3] = phigh[i] - 2;
        }
      }
      else
      {
        // For cell and face meshes always subtract one because the
        // domain is inclusive.
        for (int i=0; i<3; i++)
        {
          extents[i+0] = plow[i];
          extents[i+3] = phigh[i] - 1;
        }
      }

      // These extents are inclusive.
      dn->SetNestingForDomain(p, my_level, childPatches[p], extents);

      // debug5 << "\tdn->SetNestingForDomain("
      //        << my_level << "," << local_patch << ") <"
      //        << extents[0] << "," << extents[1] << "," << extents[2] << "> to <"
      //        << extents[3] << "," << extents[4] << "," << extents[5] << ">";

      // debug5 << "\t children patches <";

      // for (int i=0; i<childPatches[p].size(); ++i)
      // {
      //   int child_level, child_patch;
      //   GetLevelAndLocalPatchNumber(childPatches[p][i], child_level, child_patch);

      //   debug5 << "(" << child_level << "," << child_patch << "),  ";

      //   debug5 << childPatches[p][i] << ",  ";
      // }

      // debug5 << ">" << std::endl;
    }

    this->mesh_domains[meshname] =
      void_ref_ptr(dn, avtStructuredDomainNesting::Destruct);

    forceMeshReload = false;
  }

  // Register these structures with the generic database so that it knows
  // to ghost out the right cells.
  cache->CacheVoidRef("any_mesh", // key MUST be called any_mesh
                      AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                      timestate, -1, this->mesh_boundaries[meshname]);
  cache->CacheVoidRef("any_mesh", // key MUST be called any_mesh
                      AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                      timestate, -1, this->mesh_domains[meshname]);

  //VERIFY we got the mesh boundary and domain in there
  void_ref_ptr vrTmp = cache->GetVoidRef("any_mesh", // MUST be called any_mesh
                            AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                            timestate, -1);
  if (*vrTmp == NULL || *vrTmp != *this->mesh_boundaries[meshname])
    throw InvalidFilesException("uda boundary mesh not registered");

  vrTmp = cache->GetVoidRef("any_mesh", // MUST be called any_mesh
                            AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                            timestate, -1);
  if (*vrTmp == NULL || *vrTmp != *this->mesh_domains[meshname])
    throw InvalidFilesException("uda domain mesh not registered");

  // visitTimer->StopTimer(t1, "avtUintahFileFormat::GetDomainBoundariesAndNesting");
#endif
}


// ***************************************************************************
//  Method: avtUintahFileFormat::ConvertTo*
//
//  Purpose:
//      Converts a data arry to ints or longs.
//
//  Arguments:
//      Input array
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************

vtkIntArray *ConvertToInt(vtkDataArray *input)
{
  vtkIntArray *iarr = vtkIntArray::New();
  iarr->SetNumberOfComponents(input->GetNumberOfComponents());
  iarr->SetNumberOfTuples(input->GetNumberOfTuples());
  for(vtkIdType i = 0; i < input->GetNumberOfTuples(); ++i)
  {
    //iarr->SetTuple(i, input->GetTuple(i));
    double foo=i;
    iarr->SetTuple(i,&foo);
  }
  return iarr;
}

vtkLongArray *ConvertToLong(vtkDataArray *input)
{
  vtkLongArray *iarr = vtkLongArray::New();
  iarr->SetNumberOfComponents(input->GetNumberOfComponents());
  iarr->SetNumberOfTuples(input->GetNumberOfTuples());
  for(vtkIdType i = 0; i < input->GetNumberOfTuples(); ++i)
  {
    //iarr->SetTuple(i, input->GetTuple(i));
    double foo=i;
    iarr->SetTuple(i,&foo);
  }
  return iarr;
}

// ***************************************************************************
//  Method: avtUintahFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
vtkDataSet *
avtUintahFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
  // int t1 = visitTimer->StartTimer();

  //debug5<<"avtUintahFileFormat::GetMesh(timestate="<<timestate<<",domain="<<domain<<",meshname="<<meshname<<std::endl;

  ActivateTimestep(timestate);

  std::string meshName(meshname);

  int level, local_patch;
  GetLevelAndLocalPatchNumber(domain, level, local_patch);

  // Particle data
  if (meshName.find("Particle_Mesh") != std::string::npos)
  {
    size_t found = meshName.find("/");
    std::string matl = meshName.substr(found + 1);

    int matlNo = -1;
    if (matl.compare("*") != 0)
      matlNo = atoi(matl.c_str());

    // int t2 = visitTimer->StartTimer();

    // Get the particle position name typically p.x
    std::string posName((*getParticlePositionName)(archive));

    ParticleDataRaw *pd = (*getParticleData)(archive, grid, level, local_patch,
                                             posName, matlNo, timestate);

    // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetMesh getParticleData");
    // Not all patches will have particle data.
    if( pd )
    {
      // Create the vtkPoints object and copy points into it.
      vtkDoubleArray *doubleArray = vtkDoubleArray::New();
      doubleArray->SetNumberOfComponents(3);
      doubleArray->SetArray(pd->data, pd->num*pd->components, 0);

      vtkPoints *points = vtkPoints::New();
      points->SetData(doubleArray);
      doubleArray->Delete();

      // Create a vtkUnstructuredGrid to contain the point cells.
      vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
      ugrid->SetPoints(points);
      points->Delete();
      ugrid->Allocate(pd->num);

      vtkIdType onevertex;
      for(int i = 0; i < pd->num; ++i)
      {
        onevertex = i;
        ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
      }

      // Don't delete pd->data - VTK owns it now!
      delete pd;

      // Try to retrieve existing cache ref
      void_ref_ptr vrTmp =
        cache->GetVoidRef(meshname, AUXILIARY_DATA_GLOBAL_NODE_IDS,
                          timestate, domain);

      vtkDataArray *pID = NULL;

      if (*vrTmp == NULL)
      {
        // Add global node ids to facilitate point cloud usage basically
        // same as GetVar(timestate, domain, "particleID");
        int matlNo = -1;
        if (matl.compare("*") != 0)
          matlNo = atoi(matl.c_str());

        ParticleDataRaw *pd = NULL;

        //debug5<<"\t(*getParticleData)...\n";

        //todo: this returns an array of doubles. Need to return
        //expected datatype to avoid unnecessary conversion.
#if LIB_VERSION_LE(UINTAH, 2,0,0)
        if( variableExists(archive, "p.particleID") )
#endif
        {
          // int t2 = visitTimer->StartTimer();

          debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

          pd = (*getParticleData)(archive, grid, level, local_patch,
                                  "p.particleID", matlNo, timestate);

          debug5<<__FILE__<<"  "<<__FUNCTION__<<"  "<< __LINE__<< std::endl;

          // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetMesh getParticleData");
        }

        //debug5 << "got particle data: "<<pd<<"\n";
        if (pd)
        {
          vtkDoubleArray *rv = vtkDoubleArray::New();
          //vtkLongArray *rv = vtkLongArray::New();
          //debug5<<"\tSetNumberOfComponents("<<pd->components<<")...\n";
          rv->SetNumberOfComponents(pd->components);

          //debug5<<"\tSetArray...\n";
          rv->SetArray(pd->data, pd->num*pd->components, 0);

          // don't delete pd->data - VTK owns it now!
          delete pd;

          // TODO: This is the unnecesary conversion, from long
          // int->double->int, to say nothing of the implicit curtailing
          // that might occur (note also: this is a VisIt bug that uses
          // ints to store particle ids rather than long ints)
          vtkIntArray *iv = ConvertToInt(rv);
          //vtkLongArray *iv=ConvertToLong(rv);
          rv->Delete(); // this should now delete pd->data

          pID = iv;
        }

        //debug5<<"read particleID ("<<pID<<")\n";
        if (pID)
        {
          //debug5<<"adding global node ids from particleID\n";
          pID->SetName("avtGlobalNodeId");
          void_ref_ptr vr = void_ref_ptr( pID , avtVariableCache::DestructVTKObject );
          cache->CacheVoidRef( meshname, AUXILIARY_DATA_GLOBAL_NODE_IDS,
                               timestate, domain, vr);

          //make sure it worked
          void_ref_ptr vrTmp =
            cache->GetVoidRef(meshname, AUXILIARY_DATA_GLOBAL_NODE_IDS,
                              timestate, domain);

          if (*vrTmp == NULL || *vrTmp != *vr)
            throw InvalidFilesException("failed to register uda particle global node");
        }
      }

      // visitTimer->StopTimer(t1, "avtUintahFileFormat::GetMesh Particle Grid");

      return ugrid;
    }
    else
      return nullptr;
  }

  // Volume patch data
  else if (meshName.find("Patch_Mesh") != std::string::npos)
  {
    //debug5<<"Calculating vtkRectilinearGrid mesh for "<<meshName<<" mesh ("<<rgrid<<").\n";

    // Make sure we have ghosting info for this mesh
    GetDomainBoundariesAndNesting(timestate, meshname);

    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    int dims[3] = {2, 2, 2}, base[3] = {0, 0, 0};

    // Get the patch bounds
    int plow[3], phigh[3];
    patchInfo.getBounds(plow, phigh, "CC_MESH");

    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    rgrid->SetDimensions(dims);

    // Set the coordinates of the grid points in each direction.
    for (int c=0; c<3; c++)
    {
      vtkFloatArray *coords = vtkFloatArray::New();
      coords->SetNumberOfTuples(dims[c]);
      float *array = (float *) coords->GetVoidPointer(0);

      for (int i=0; i<dims[c]; i++)
      {
        array[0] = levelInfo.anchor[c] +  plow[c] * levelInfo.spacing[c];
        array[1] = levelInfo.anchor[c] + phigh[c] * levelInfo.spacing[c];
      }

      switch(c)
      {
      case 0:
        rgrid->SetXCoordinates(coords);
        break;
      case 1:
        rgrid->SetYCoordinates(coords);
        break;
      case 2:
        rgrid->SetZCoordinates(coords);
        break;
      }

      coords->Delete();
    }

    // visitTimer->StopTimer(t1, "avtUintahFileFormat::GetMesh Volume Patch Grid");

    return rgrid;
  }

  // Volume grid data
  else
  {
    //debug5<<"Calculating vtkRectilinearGrid mesh for "<<meshName<<" mesh ("<<rgrid<<").\n";

    // Make sure we have ghosting info for this mesh
    GetDomainBoundariesAndNesting(timestate, meshname);

    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    int dims[3];

    // Get the patch bounds
    int plow[3], phigh[3];
    patchInfo.getBounds(plow, phigh, meshName);

    // For node based meshes add one if there is a neighbor patch.
    if( meshName.find("NC_") == 0 )
    {
      int nlow[3], nhigh[3];
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

      for (int i=0; i<3; i++)
      {
        phigh[i] += nhigh[i];
        dims[i] = phigh[i] - plow[i];
      }
    }
    else
    {
      // For cell and face meshes always add one.
      for (int i=0; i<3; i++)
      {
        dims[i] = phigh[i] - plow[i] + 1;
      }
    }

    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    rgrid->SetDimensions(dims);

    // These are needed to offset grid points in order to preserve
    // face centered locations on node-centered domain.
    bool sfck[3] = { meshName.find("SFCX") != std::string::npos,
                     meshName.find("SFCY") != std::string::npos,
                     meshName.find("SFCZ") != std::string::npos };

    int nlow[3], nhigh[3];

    if( sfck[0] || sfck[1] || sfck[2] )
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

    // Set the coordinates of the grid points in each direction.
    for (int c=0; c<3; c++)
    {
      vtkFloatArray *coords = vtkFloatArray::New();
      coords->SetNumberOfTuples(dims[c]);
      float *array = (float *) coords->GetVoidPointer(0);

      for (int i=0; i<dims[c]; i++)
      {
        // Face centered data gets shifted towards -inf by half a cell.
        // Boundary patches are special shifted to preserve global domain.
        // Internal patches are always just shifted.
        float face_offset;

        if (sfck[c])
        {
          if (i==0)
          {
            // No neighbor, so the patch is on low boundary
            if (nlow[c] == 0)
              face_offset = 0.0;
            // Patch boundary is internal to the domain
            else
              face_offset = -0.5;
          }
          else if (i == dims[c]-1)
          {
            // No neighbor, so the patch is on high boundary
            if (nhigh[c] == 0)
            {
              // Periodic means one less value in the face-centered direction
              if (levelInfo.periodic[c])
                face_offset = 0.0;
              else
                face_offset = -1;
            }
            // Patch boundary is internal to the domain
            else
              face_offset = -0.5;
          }
          else
            face_offset = -0.5;
        }
        else
          face_offset = 0;

        array[i] = levelInfo.anchor[c] +
          (i + plow[c] + face_offset) * levelInfo.spacing[c];
      }

      switch(c)
      {
      case 0:
        rgrid->SetXCoordinates(coords);
        break;
      case 1:
        rgrid->SetYCoordinates(coords);
        break;
      case 2:
        rgrid->SetZCoordinates(coords);
        break;
      }

      coords->Delete();
    }

    // visitTimer->StopTimer(t1, "avtUintahFileFormat::GetMesh Volume Grid");

    return rgrid;
  }

  return NULL;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
vtkDataArray *
avtUintahFileFormat::GetVar(int timestate, int domain, const char *varname)
{
  // int t1 = visitTimer->StartTimer();

  ActivateTimestep(timestate);

  bool isParticleVar = false;
  bool isInternalVar = false;

  // Get the var name sans the material. If a patch variable then the
  // var name will be "Patch".
  std::string varName(varname);
  size_t found = varName.find("/");
  std::string matl = varName.substr(found + 1);
  varName = varName.substr(0, found);

  // Get the varType except for patch based data which does not come
  // from the data warehouse but instead from internal structures.
  std::string varType;

  // Patch based data ids and bounds
  if( strncmp(varname, "Patch/Id", 8) == 0 ||
      strncmp(varname, "Patch/ProcId", 12) == 0 ||
      strncmp(varname, "Patch/Bounds/Low",  16) == 0 ||
      strncmp(varname, "Patch/Bounds/High", 17) == 0 )
  {
    isInternalVar = true;
    varType = "Patch_Mesh";
  }
  // Patch based node ids
  else if( strncmp(varname, "Patch/Nodes", 11) == 0 )
  {
    isInternalVar = true;
    varType = "NC_Mesh";
  }
  // Patch node and ranks.
  else if( strncmp(varname, "Nodes/", 6) == 0 ||
           strncmp(varname, "Ranks/", 6) == 0 )
  {
    isInternalVar = true;

    if(strncmp(&(varname[6]), "Particle_Mesh", 13) == 0)
    {
      isParticleVar = true;

      varType = "Particle_Mesh";

      // Get the material sans "Particle_Mesh".
      varName = std::string(varname);
      found = varName.find_last_of("/");
      matl = varName.substr(found + 1);
    }
    else
    {
      // Get the mesh type.
      varName = std::string(varname);
      found = varName.find_last_of("/");
      varType = varName.substr(found + 1);
    }
  }
  // Grid variables
  else
  {
    // For PerPatch data remove the Patch/ prefix and get the var
    // name and set the material to zero.
    if( varName == "Patch" )
    {
      // Get the var name and material sans "Patch/".
      varName = std::string(varname);
      found = varName.find("/");
      varName = varName.substr(found + 1);

      matl = "0";
    }

    // Get the var type and check to see if it is a particle variable.
    for (int k=0; k<(int)stepInfo->varInfo.size(); k++)
    {
      if (stepInfo->varInfo[k].name == varName)
      {
        varType = stepInfo->varInfo[k].type;

        // Check for a particle variable
        if (stepInfo->varInfo[k].type.find("ParticleVariable") !=
            std::string::npos)
        {
          isParticleVar = true;
        }

        break;
      }
    }
  }

  vtkDoubleArray *rv = nullptr;

  int level, local_patch;
  GetLevelAndLocalPatchNumber(domain, level, local_patch);

  // particle data
  if (isParticleVar)
  {
    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    int matlNo = -1;
    if (matl.compare("*") != 0)
      matlNo = atoi(matl.c_str());

    ParticleDataRaw *pd = NULL;

    if( isInternalVar )
    {
      // Patch node and ranks.
      if( strncmp(varname, "Nodes/", 6) == 0 ||
          strncmp(varname, "Ranks/", 6) == 0 )
      {
        unsigned int numParticles =
          getNumberParticles(archive, grid, level, local_patch, matlNo, timestate);
        // Some patches may not have particles.
        if( numParticles )
        {
          pd = new ParticleDataRaw;
          pd->num = numParticles;
          pd->components = 1;
          pd->data = new double[pd->num * pd->components];

          // Patch processor rank
#if LIB_VERSION_LE(UINTAH, 2,5,1)
          double value = patchInfo.getProcId();
#elif LIB_VERSION_LE(UINTAH, 2,2,0)
          double value = patchInfo.getProcRankId();
#else
          double value = patchInfo.getProcId();
#endif
          for (int i=0; i<pd->num*pd->components; ++i)
            pd->data[i] = value;
        }
      }
    }
    else
    {
#ifdef SERIALIZED_READS
      int numProcs, rank;
      int msg = 128, tag = 256;
      MPI_Status status;

      MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
      MPI_Comm_rank(VISIT_MPI_COMM, &rank);

      int totalPatches = 0;
      for (int i = 0; i < stepInfo->levelInfo.size(); i++)
        totalPatches += stepInfo->levelInfo[i].patchInfo.size();

      // Calculate which process we should wait for a message from
      // if we're processing doiman 0 don't wait for anyone else
      int prev = (rank+numProcs-1)%numProcs;
      int next = (rank+1)%numProcs;

      // Domain 0 will always reads right away
      if (domain == 0)
        prev = -1;
      //debug5 << "Proc: " << rank << " sent to GetVar" << endl;

      // Wait for the previous read to finish
      if (prev>=0)
        MPI_Recv(&msg, 1, MPI_INT, prev, tag, VISIT_MPI_COMM, &status);

      // int t2 = visitTimer->StartTimer();

      pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);

      // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getParticleData");

      // Let the next read go
      if (next>=0)
        MPI_Send(&msg, 1, MPI_INT, next, tag, VISIT_MPI_COMM);

#else
      // int t2 = visitTimer->StartTimer();

      pd = (*getParticleData)(archive, grid, level, local_patch, varName, matlNo, timestate);

      // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getParticleData");
#endif

      if( pd )
        CheckNaNs(pd->data, pd->num*pd->components, varname, level, local_patch);
    }

    // Some patches may not have particles.
    if( pd )
    {
      rv = vtkDoubleArray::New();
      rv->SetNumberOfComponents(pd->components);
      rv->SetArray(pd->data, pd->num*pd->components, 0);

      // don't delete pd->data - vtk owns it now!
      delete pd;
    }
  }

  // Volume data
  else
  {
    LevelInfo &levelInfo = stepInfo->levelInfo[level];
    PatchInfo &patchInfo = levelInfo.patchInfo[local_patch];

    // The region we're going to ask uintah for (from plow to phigh-1)
    int plow[3], phigh[3];
    patchInfo.getBounds(plow, phigh, varType);

    // For node based meshes add one if there is a neighbor patch.
    bool nodeCentered = (varType.find("NC") != std::string::npos);

    if( nodeCentered )
    {
      int nlow[3], nhigh[3];
      patchInfo.getBounds(nlow, nhigh, "NEIGHBORS");

      for (int i=0; i<3; i++)
        phigh[i] += nhigh[i];
    }

#if LIB_VERSION_LE(UINTAH, 2,5,1)
    double rank = patchInfo.getProcId();
#elif LIB_VERSION_LE(UINTAH, 2,2,0)
    double rank = patchInfo.getProcRankId();
#else
    double rank = patchInfo.getProcId();
#endif

    GridDataRaw *gd = NULL;

    // The data for these variables does not come from the data
    // warehouse but instead from internal structures.
    if( isInternalVar )
    {
      gd = new GridDataRaw;

      // Patch ids.
      if (strcmp(varname, "Patch/Id") == 0) {
        gd->num = 1;
        gd->components = 1;
        gd->data = new double[gd->num * gd->components];

#if LIB_VERSION_LE(UINTAH, 2,0,0)
        gd->data[0] = patchInfo.getPatchId();
#endif
      }
      // Patch processor ids.
      else if (strcmp(varname, "Patch/ProcId") == 0 ) {
        gd->num = 1;
        gd->components = 1;
        gd->data = new double[gd->num * gd->components];
        gd->data[0] = rank;
      }
      // Patch based bounds.
      else if( strncmp(varname, "Patch/Bounds/Low",  16) == 0 ||
               strncmp(varname, "Patch/Bounds/High", 17) == 0) {
        gd->num = 1;
        gd->components = 3; // Bounds are vectors
        gd->data = new double[gd->num * gd->components];

        // Get the bounds for this mesh as a variable (not for the grid).
        std::string meshname = std::string(varname);
        found = meshname.find_last_of("/");
        meshname = meshname.substr(found + 1);

        patchInfo.getBounds(plow, phigh, meshname);

        int *value;
        if (strncmp(varname, "Patch/Bounds/Low", 16) == 0 )
          value = &plow[0];
        else //if( strncmp(varname, "Patch/Bounds/High", 17) == 0)
          value = &phigh[0];

        for (int c=0; c<3; c++)
          gd->data[c] = value[c];
      }
      // Patch based node ids
      else if( strncmp(varname, "Patch/Nodes",  11) == 0 ) {
        // Using the node mesh on the simulation
        gd->num = ((phigh[0] - plow[0]) *
                   (phigh[1] - plow[1]) *
                   (phigh[2] - plow[2]));
        // Bounds are vectors
        gd->components = 3;
        gd->data = new double[gd->num * gd->components];

        int cc = 0;

        for( int k=plow[2]; k<phigh[2]; ++k )
        {
          for( int j=plow[1]; j<phigh[1]; ++j )
          {
            for( int i=plow[0]; i<phigh[0]; ++i )
            {
              gd->data[cc++] = i;
              gd->data[cc++] = j;
              gd->data[cc++] = k;
            }
          }
        }
      }
      // Patch based data node and ranks.
      else if( strncmp(varname, "Nodes/", 6) == 0 ||
               strncmp(varname, "Ranks/", 6) == 0 ) {

        // Using the patch mesh
        if (strcmp(&(varname[6]), "Patch_Mesh") == 0 )
          gd->num = 1;
        // Using the cell mesh.
        else /* if (strcmp(&(varname[6]), "CC_Mesh") == 0 ||
                    strcmp(&(varname[6]), "NC_Mesh") == 0 ||
                    strcmp(&(varname[6]), "SFCX_Mesh") == 0 ||
                    strcmp(&(varname[6]), "SFCY_Mesh") == 0 ||
                    strcmp(&(varname[6]), "SFCZ_Mesh") == 0 ) */
          gd->num = ((phigh[0] - plow[0]) *
                     (phigh[1] - plow[1]) *
                     (phigh[2] - plow[2]));

        gd->components = 1;
        gd->data = new double[gd->num * gd->components];

        for (int i=0; i<gd->num; i++)
          gd->data[i] = rank;
      }
    }

    // Patch grid data from the warehouse
    else
    {
      int matlNo = -1;
      if (matl.compare("All") != 0)
        matlNo = atoi(matl.c_str());

#ifdef SERIALIZED_READS
      int numProcs, rank;
      int msg = 128, tag = 256;
      MPI_Status status;

      MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
      MPI_Comm_rank(VISIT_MPI_COMM, &rank);

      int totalPatches = 0;
      for (int i = 0; i < stepInfo->levelInfo.size(); i++)
        totalPatches += stepInfo->levelInfo[i].patchInfo.size();

      // calculate which process we should wait for a message from
      // if we're processing doiman 0 don't wait for anyone else
      int prev = (rank+numProcs-1)%numProcs;
      int next = (rank+1)%numProcs;

      // domain 0 always reads right away
      if (domain==0)
        prev = -1;
      //debug5 << "Proc: " << rank << " sent to GetVar" << endl;

      // wait for previous read to finish
      if (prev>=0)
        MPI_Recv(&msg, 1, MPI_INT, prev, tag, VISIT_MPI_COMM, &status);

      // int t2 = visitTimer->StartTimer();

      gd = (*getGridData)(archive, grid, level, local_patch,
                          varName, matlNo, timestate,
                          plow, phigh,
                          (nodeCentered ? 0 : loadExtraGeometry));

      // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getGridData");

      // let the next read go
      if (next>=0)
        MPI_Send(&msg, 1, MPI_INT, next, tag, VISIT_MPI_COMM);
#else
      // int t2 = visitTimer->StartTimer();

      gd = (*getGridData)(archive, grid, level, local_patch,
                          varName, matlNo, timestate,
                          plow, phigh,
                          (nodeCentered ? 0 : loadExtraGeometry));

      // visitTimer->StopTimer(t2, "avtUintahFileFormat::GetVar getGridData");
#endif
    }

    if( gd  )
    {
      CheckNaNs(gd->data, gd->num*gd->components, varname, level, local_patch);

      rv = vtkDoubleArray::New();
      rv->SetNumberOfComponents(gd->components);
      rv->SetArray(gd->data, gd->num*gd->components, 0);

      // Don't delete gd->data - VTK owns it now!
      delete gd;
    }
    else
    {
      std::stringstream msg;
      msg << "Visit - "
          << "Uintah variable \"" << varname << "\"  "
          << "could not be processed.";

      avtCallback::IssueWarning(msg.str().c_str());
    }
  }

  // visitTimer->StopTimer(t1, "avtUintahFileFormat::GetVar");

  return rv;
}


// ****************************************************************************
//  Method: avtUintahFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: sshankar -- generated by xml2avt
//  Creation:   Tue May 13 19:02:26 PST 2008
//
// ****************************************************************************
vtkDataArray *
avtUintahFileFormat::GetVectorVar(int timestate, int domain,const char *varname)
{
  // Handle vector variables exactly the same way as scalar variables
  return GetVar(timestate, domain, varname);
}


// ****************************************************************************
//  Method: avtBoxlib2DFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data specified.
//
//  Arguments:
//      var        The variable of interest.
//      dom        The domain of interest.
//      type       The type of auxiliary data.
//      <unnamed>  The arguments for that type -- not used.
//      df         Destructor function.
//
//  Returns:    The auxiliary data.
//
//  Programmer: Hank Childs
//  Creation:   January 22, 2006
//
// ****************************************************************************
void *
avtUintahFileFormat::GetAuxiliaryData(const char *var, int dom,
                                      const char * type, void *,
                                      DestructorFunction &df)
{
  return NULL;
}

// ****************************************************************************
//  Method: avtUintahFileFormat::CheckNaNs
//
//  Purpose:
//      Check for and warns about NaN values in the data.
//
//  Arguments:
//      num        data size
//      data       data
//      level      level that contains this patch
//      patch      patch that contains these cells
//
//  Returns:    none
//
//  Programmer: cchriste
//  Creation:   06.02.2012
//
// ****************************************************************************
void
avtUintahFileFormat::CheckNaNs(double *data, const int num,
                               const char* varname,
                               const int level, const int patch)
{
  // Replace nan's with a large negative number
  std::vector<int> nanCells;

  for (int i=0; i<num; i++)
  {
    if (std::isnan(data[i]))
    {
      data[i] = NAN_REPLACE_VAL;
      nanCells.push_back(i);
    }
  }

  if (!nanCells.empty())
  {
    std::stringstream sstr;
    sstr << "NaNs exist for variable " << varname
         << " in patch " << patch << " of level " << level
         << " and " << nanCells.size() << "/" << num
         << " cells have been replaced by the value "
         <<  NAN_REPLACE_VAL << ".";

    // if ((int)nanCells.size()>40)
    // {
    //   sstr<<"\nFirst 20: ";
    //   for (int i=0;i<(int)nanCells.size() && i<20;i++)
    //     sstr<<nanCells[i]<<",";
    //   sstr<<"\nLast 20: ";
    //   for (int i=(int)nanCells.size()-21;i<(int)nanCells.size();i++)
    //     sstr<<nanCells[i]<<",";
    // }
    // else
    // {
    //   for (int i=0;i<(int)nanCells.size();i++)
    //     sstr<<nanCells[i]<<((int)nanCells.size()!=(i+1)?",":".");
    // }

    avtCallback::IssueWarning(sstr.str().c_str());
  }
}


// ****************************************************************************
//  Method: avtUintahFileFormat::addRectilinearMesh
//
//  Purpose:
//      Add a rectilinear mesh
//
//  Arguments:
//      md        Database meta data
//
//  Returns:    none
//
//  Programmer: Allen R. Sanderson
//  Creation:   12 Dec. 2018
//
// ****************************************************************************
void
avtUintahFileFormat::addRectilinearMesh( avtDatabaseMetaData *md,
                                         std::set<std::string> &meshes_added,
                                         std::string meshName,
                                         TimeStepInfo* stepInfo )
{
  if (meshes_added.find(meshName) == meshes_added.end())
  {
    int numLevels = stepInfo->levelInfo.size();

    int totalPatches = 0;
    for (int i = 0; i < numLevels; i++)
      totalPatches += stepInfo->levelInfo[i].patchInfo.size();

    // Used for the domain partitioning
    std::vector<std::string> pieceNames(totalPatches);

    for (int i = 0; i < totalPatches; i++)
      {
        char tmpName[64];
        int level, local_patch;

        GetLevelAndLocalPatchNumber(i, level, local_patch);
        sprintf(tmpName,"level%d, patch%d", level, local_patch);

        pieceNames[i] = tmpName;
      }

    // Compute the bounding box of the mesh from the grid indices of
    // level 0.
    LevelInfo &levelInfo = stepInfo->levelInfo[0];

    // This can be done once for everything because the spatial range is
    // the same for all meshes.
    double box_min[3], box_max[3];
    levelInfo.getExtents(box_min, box_max);

    avtMeshMetaData *mesh = new avtMeshMetaData;

    mesh->name = meshName;
    mesh->meshType = AVT_AMR_MESH;
    mesh->topologicalDimension = 3;
    mesh->spatialDimension = 3;

    mesh->numBlocks = totalPatches;
    mesh->blockTitle = "patches";
    mesh->blockPieceName = "patch";
    mesh->numGroups = numLevels;
    mesh->groupTitle = "levels";
    mesh->groupPieceName = "level";
    mesh->blockNames = pieceNames;
    mesh->containsExteriorBoundaryGhosts = false;

    mesh->hasSpatialExtents = true;
    mesh->minSpatialExtents[0] = box_min[0];
    mesh->maxSpatialExtents[0] = box_max[0];
    mesh->minSpatialExtents[1] = box_min[1];
    mesh->maxSpatialExtents[1] = box_max[1];
    mesh->minSpatialExtents[2] = box_min[2];
    mesh->maxSpatialExtents[2] = box_max[2];

    int low[3], high[3];
    levelInfo.getBounds(low, high, meshName);
    int logical[3];
    for (int i=0; i<3; ++i)
      logical[i] = high[i]-low[i];

    mesh->hasLogicalBounds = true;
    mesh->logicalBounds[0] = logical[0];
    mesh->logicalBounds[1] = logical[1];
    mesh->logicalBounds[2] = logical[2];

    md->Add(mesh);
    meshes_added.insert(meshName);

    addMeshNodeRankSIL( md, meshName );

    meshes_added.insert(meshName);
  }
}

// ****************************************************************************
//  Method: avtUintahFileFormat::addParticleMesh
//
//  Purpose:
//      Add a particle mesh
//
//  Arguments:
//      md        Database meta data
//
//  Returns:    none
//
//  Programmer: Allen R. Sanderson
//  Creation:   12 Dec. 2018
//
// ****************************************************************************
void
avtUintahFileFormat::addParticleMesh( avtDatabaseMetaData *md,
                                      std::set<std::string> &meshes_added,
                                      std::string meshName,
                                      TimeStepInfo* stepInfo )
{
  if (meshes_added.find(meshName) == meshes_added.end())
  {
    int numLevels = stepInfo->levelInfo.size();

    int totalPatches = 0;
    for (int i = 0; i < numLevels; i++)
      totalPatches += stepInfo->levelInfo[i].patchInfo.size();

    // Used for the domain partitioning
    std::vector<std::string> pieceNames(totalPatches);

    for (int i = 0; i < totalPatches; i++)
      {
        char tmpName[64];
        int level, local_patch;

        GetLevelAndLocalPatchNumber(i, level, local_patch);
        sprintf(tmpName,"level%d, patch%d", level, local_patch);

        pieceNames[i] = tmpName;
      }

    // Compute the bounding box of the mesh from the grid indices of
    // level 0.
    LevelInfo &levelInfo = stepInfo->levelInfo[0];

    // This can be done once for everything because the spatial range is
    // the same for all meshes.
    double box_min[3], box_max[3];
    levelInfo.getExtents(box_min, box_max);

    avtMeshMetaData *mesh = new avtMeshMetaData;

    mesh->name = meshName;
    mesh->meshType = AVT_POINT_MESH;
    mesh->topologicalDimension = 0;
    mesh->spatialDimension = 3;

    mesh->numBlocks = totalPatches;
    mesh->blockTitle = "patches";
    mesh->blockPieceName = "patch";
    mesh->numGroups = numLevels;
    mesh->groupTitle = "levels";
    mesh->groupPieceName = "level";
    mesh->blockNames = pieceNames;

    mesh->hasSpatialExtents = true;
    mesh->minSpatialExtents[0] = box_min[0];
    mesh->maxSpatialExtents[0] = box_max[0];
    mesh->minSpatialExtents[1] = box_min[1];
    mesh->maxSpatialExtents[1] = box_max[1];
    mesh->minSpatialExtents[2] = box_min[2];
    mesh->maxSpatialExtents[2] = box_max[2];

    int low[3], high[3];
    levelInfo.getBounds(low, high, meshName);
    int logical[3];
    for (int i=0; i<3; ++i)
      logical[i] = high[i]-low[i];

    // mesh->hasLogicalBounds = true;
    // mesh->logicalBounds[0] = logical[0];
    // mesh->logicalBounds[1] = logical[1];
    // mesh->logicalBounds[2] = logical[2];

    md->Add(mesh);
    meshes_added.insert(meshName);

    addMeshNodeRankSIL( md, meshName );

    meshes_added.insert(meshName);
  }
}

// ****************************************************************************
//  Method: avtUintahFileFormat::addMeshNodeRankSIL
//
//  Purpose:
//      Adds a mesh SIL for nodes and ranks
//
//  Arguments:
//      md        Database meta data
//      mesh      Mesh for the SIL
//
//  Returns:    none
//
//  Programmer: Allen R. Sanderson
//  Creation:   12 Dec. 2018
//
// ****************************************************************************
void
avtUintahFileFormat::addMeshNodeRankSIL( avtDatabaseMetaData *md,
                                         std::string mesh )
{
  // Add a SIL for subsettng via the nodes and ranks.
  const unsigned int nRanksPerGroup = ceil((float) nProcs / 10.0);
  const unsigned int nRanks = nProcs;
  const unsigned int nNodes = nRanksPerGroup > 5 ?
    ceil((float)nRanks/(float)nRanksPerGroup) : 1;

  int rank_enum_id[nRanks];
  int node_enum_id[nNodes];

  std::string enum_name = std::string("Ranks/") + mesh;

  avtScalarMetaData *smd =
    new avtScalarMetaData(enum_name, mesh, AVT_ZONECENT );

  smd->SetEnumerationType(avtScalarMetaData::ByValue);
  smd->hideFromGUI = true;

  for( unsigned int i=0; i<nRanks; ++i ) {
    char msg[128];
    sprintf( msg, "Rank_%04d", i );
    rank_enum_id[i] = smd->AddEnumNameValue( msg, i);
  }

  if( nNodes > 1 ) {

    for( unsigned int i=0, j=1; i<nNodes; ++i, ++j ) {
      char msg[128];
      sprintf( msg, "Ranks_%04d_%04d",
               i*nRanksPerGroup,
               std::min(j*nRanksPerGroup, nProcs)-1 );
      node_enum_id[i] = smd->AddEnumNameValue( msg, nRanks+i);
    }

    for( unsigned int i=0; i<nRanks; ++i ) {
      smd->AddEnumGraphEdge(node_enum_id[i/nRanksPerGroup], rank_enum_id[i], "Ranks" );
    }
  }

  md->Add(smd);
}

// ****************************************************************************
//  Method: avtUintahFileFormat::addMeshVariable
//
//  Purpose:
//      Adds a mesh variable
//
//  Arguments:
//      md        Database meta data
//      mesh      Mesh for the SIL
//
//  Returns:    none
//
//  Programmer: Allen R. Sanderson
//  Creation:   12 Dec. 2018
//
// ****************************************************************************
void
avtUintahFileFormat::addMeshVariable( avtDatabaseMetaData *md,
                                      std::set<std::string> &mesh_vars_added,
                                      std::string varName, std::string varType,
                                      std::string meshName, avtCentering cent )
{
  // Make sure a variable is added only once.
  if (mesh_vars_added.find(meshName+varName) == mesh_vars_added.end())
  {
    mesh_vars_added.insert(meshName+varName);

    // 3 -> point/vector dimension
    if (varType.find("Point") != std::string::npos ||
        varType.find("Vector") != std::string::npos ||
        varType.find("IntVector") != std::string::npos)
    {
      AddVectorVarToMetaData(md, varName, meshName, cent, 3);
    }
    // 9 -> tensor
    else if (varType.find("Matrix3") != std::string::npos)
    {
      AddTensorVarToMetaData(md, varName, meshName, cent, 9);
    }
    // 7 -> vector
    else if (varType.find("Stencil7") != std::string::npos)
    {
      AddVectorVarToMetaData(md, varName, meshName, cent, 7);
    }
    // 4 -> vector
    else if (varType.find("Stencil4") != std::string::npos)
    {
      AddVectorVarToMetaData(md, varName, meshName, cent, 4);
    }
    // 1 -> scalar
    else
    {
      AddScalarVarToMetaData(md, varName, meshName, cent);
    }
    // else
    // {
    //   std::stringstream msg;
    //   msg << "Visit - "
    //       << "Uintah variable \"" << varName << "\"  "
    //       << "has an unknown variable type \""
    //       << varType << "\"";

    //   avtCallback::IssueWarning(msg.str().c_str());

    //   return;
    // }
  }
}
